
################################################
# Media Influence and Spatial Voting: The Role of Perceived Party Positions
# Political Behavior
# Lucas Paulo da Silva
# RScript for performing main analysis
# 28/02/2025
################################################


##### Set up ######

# R version 4.4.0 (2024-04-24 ucrt)

# load
library(haven)
library(estimatr)
library(dplyr)
library(ggplot2)
library(texreg)
library(ggtext)
library(plotrix)
library(dotwhisker)
library(tidyr)

require(rio)
require(foreign)
require(car)
require(plyr)
require(plm)
require(lme4)
library(multiwayvcov)
require(lmtest)
require(stargazer)
require(xtable)
require(psych)
require(grDevices)
require(plyr)
require(reshape2)
require(pBrackets)
require(MatchIt)
require(sandwich)

# clear environment
rm(list=ls())

# set working directory
setwd("C:/your/working/directory/")

# data
bsa <-read_dta("FB_BSA.DTA")
bsa <- bsa[bsa$North==1,] #remove non-northern observations
bsa <- bsa[!is.na(bsa$year),] #drop observations without a year

# reset working directory if you want a different folder for saved plots below
#setwd("C:/your/working/directory/")

# transform class variable
bsa$lmh_l <- as.factor(dplyr::recode(as.factor(bsa$lmh), 
                            "0" = 1,"1"=0,"2"=0))

# make liverpool, treatment, and constituency variables into factors
bsa$liverpool <- as.factor(bsa$liverpool)
bsa$treat <- as.factor(bsa$treat)
bsa$constituency97_id <- as.factor(bsa$constituency97_id)

# remove years after Sun endorsed Labour
bsa <- bsa[bsa$year <= 1996,]





##### Demonstrating change in Sun readership ######

# overall
mean(bsa$sun[bsa$treat==0 & bsa$liverpool==1],na.rm=T)*100 #Liverpool, before
mean(bsa$sun[bsa$treat==1 & bsa$liverpool==1],na.rm=T)*100 #Liverpool, after

mean(bsa$sun[bsa$treat==0 & bsa$liverpool==0],na.rm=T)*100 #Northern England (NE), before
mean(bsa$sun[bsa$treat==1 & bsa$liverpool==0],na.rm=T)*100 #NE, after

# by class in Liverpool
mean(bsa$sun[bsa$treat==0 & bsa$liverpool==1 & bsa$lmh_l==1],na.rm=T)*100 #LWC
mean(bsa$sun[bsa$treat==1 & bsa$liverpool==1 & bsa$lmh_l==1],na.rm=T)*100 #LWC

mean(bsa$sun[bsa$treat==0 & bsa$liverpool==1 & bsa$lmh_l==0],na.rm=T)*100 #MHC
mean(bsa$sun[bsa$treat==1 & bsa$liverpool==1 & bsa$lmh_l==0],na.rm=T)*100 #MHC





##### Descriptive data figure ######

### general

# only use data before treatment
bsa_preboycott <- bsa[bsa$treat==1,]

# create function for data summary
data_summary <- function(data, varname, groupnames){
  require(plyr)
  require(plotrix)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE),
      se = std.error(x[[col]], na.rm=TRUE))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
  return(data_sum)
}



### panel A: perceived party position (PPP) by location

# create function for line plot ticks
number_ticks <- function(n) {function(limits) pretty(limits, n)}

# create dataframe
desc_df_labextr <- data_summary(bsa_preboycott[!(bsa_preboycott$labextr_p=="NA" | bsa_preboycott$liverpool=="NA"),], 
                                varname="labextr_p", groupnames="liverpool")

# drop NaN observation
desc_df_labextr <- desc_df_labextr[1:2,]

# create plot
ggplot(desc_df_labextr, aes(x=liverpool, y=labextr_p, fill=liverpool)) +
  geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=labextr_p-se*1.96, ymax=labextr_p+se*1.96), width=.04, size=1) +
  scale_fill_manual(values = c("0"="#6baed6", "1"="#08519c"),
                    labels = c("Northern England","Liverpool")) +
  scale_y_continuous(limits = c(0, 1),breaks=number_ticks(5)) +
  scale_x_discrete(labels=c("0" = "Northern England", "1" = "Liverpool")) +
  labs(x="Region",
       y="Proportion See Labour as Extreme",
       fill="Key") +
  geom_text(aes(label=round(labextr_p,2)), vjust=3, color="white",
            position = position_dodge(0.9), size=5, fontface = "bold")+
  theme(
    axis.title.y = element_text(size = 16),
    axis.text.y = element_text(size = 14),
    axis.title.x = element_text(size = 16),
    axis.text.x = element_text(size = 14, vjust=0.7),
    title = element_text(size = 16),
    legend.position = "none")

# save
ggsave("descriptive_ppp.png", width = 5, height = 4, units = "in", dpi = 300)



### panel B: ideological position by location

# create dataframe
desc_df_unions <- data_summary(bsa_preboycott[!(bsa_preboycott$unions_lr=="NA" | bsa_preboycott$liverpool=="NA"),], 
                               varname="unions_lr", groupnames="liverpool")

# drop NaN observation
desc_df_unions <- desc_df_unions[1:2,]

# create bar plot
ggplot(desc_df_unions, aes(x=liverpool, y=unions_lr, fill=liverpool)) +
  geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=unions_lr-se*1.96, ymax=unions_lr+se*1.96), width=.04, size=1) +
  scale_fill_manual(values = c("0"="#6baed6", "1"="#08519c"),
                    labels = c("Northern England","Liverpool")) +
  scale_y_continuous(limits = c(0, 1),breaks=number_ticks(5)) +
  scale_x_discrete(labels=c("0" = "Northern England", "1" = "Liverpool")) +
  labs(x="Region",
       y="Proportion Against Union Power",
       fill="Key") +
  geom_text(aes(label=round(unions_lr,2)), vjust=3, color="white",
            position = position_dodge(0.9), size=5, fontface = "bold")+
  theme(
    axis.title.y = element_text(size = 16),
    axis.text.y = element_text(size = 14),
    axis.title.x = element_text(size = 16),
    axis.text.x = element_text(size = 14, vjust=0.7),
    title = element_text(size = 16),
    legend.position = "none")

# save
ggsave("descriptive_ideology.png", width = 5, height = 4, units = "in", dpi = 300)



### panel C: party support by location

# create dataframe
desc_df_labsup <- data_summary(bsa_preboycott[!(bsa_preboycott$labour=="NA" | bsa_preboycott$liverpool=="NA"),], varname="labour", 
                               groupnames="liverpool")

# create bar plot
ggplot(desc_df_labsup, aes(x=liverpool, y=labour, fill=liverpool)) +
  geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=labour-se*1.96, ymax=labour+se*1.96), width=.04, size=1) +
  scale_fill_manual(values = c("0"="#6baed6", "1"="#08519c"),
                    labels = c("Northern England","Liverpool")) +
  scale_y_continuous(limits = c(0, 1),breaks=number_ticks(5)) +
  scale_x_discrete(labels=c("0" = "Northern England", "1" = "Liverpool")) +
  labs(x="Region",
       y="Proportion that Support Labour",
       fill="Key") +
  geom_text(aes(label=round(labour,2)), vjust=3, color="white",
            position = position_dodge(0.9), size=5, fontface = "bold")+
  theme(
    axis.title.y = element_text(size = 16),
    axis.text.y = element_text(size = 14),
    axis.title.x = element_text(size = 16),
    axis.text.x = element_text(size = 14, vjust=0.7),
    title = element_text(size = 16),
    legend.position = "none")

# save
ggsave("descriptive_support.png", width = 5, height = 4, units = "in", dpi = 300)





##### Comparing LWC to others ######

# percentage of LWC
(nrow(bsa[bsa$lmh_l == 1,]) / nrow(bsa)) * 100

# imp_age, political attitudes, education, voteleave, employment

# age
mean(bsa$imp_age[bsa$lmh_l == 0], na.rm = TRUE) #non-LWC
mean(bsa$imp_age[bsa$lmh_l == 1], na.rm = TRUE) #non-LWC

# female
mhc_female <- mean(bsa$female[bsa$lmh_l == 0], na.rm = TRUE) #non-LWC
print(mhc_female*100)
lwc_female <- mean(bsa$female[bsa$lmh_l == 1], na.rm = TRUE) #LWC
print(lwc_female*100)
(lwc_female - mhc_female) *100

# ethnicity
mhc_white <- mean(bsa$white[bsa$lmh_l == 0], na.rm = TRUE) #non-LWC
print(mhc_white*100)
lwc_white<- mean(bsa$white[bsa$lmh_l == 1], na.rm = TRUE) #LWC
print(lwc_white*100)
(lwc_white - mhc_white) *100

# own home
mhc_ownhome <- mean(bsa$ownhome[bsa$lmh_l == 0], na.rm = TRUE) #non-LWC
print(mhc_ownhome*100)
lwc_ownhome <- mean(bsa$ownhome[bsa$lmh_l == 1], na.rm = TRUE) #LWC
print(lwc_ownhome*100)
(lwc_ownhome - mhc_ownhome) *100

# education
bsa <- bsa %>% mutate(education_bin = case_when(
  education <= 2 ~ "Low", #ended at 16 or younger
  education >= 3 ~ "High")) #ended at 17 or older

edu_table <- prop.table(table(bsa$lmh_l, bsa$education_bin), margin=1)
print(edu_table*100)
(edu_table[1] - edu_table[2]) * 100

# religion
rel_table <- prop.table(table(bsa$lmh_l, bsa$religion), margin=1)
print(rel_table*100)
(rel_table[1] - rel_table[2]) * 100





##### Placebo tests ######

### data
prebsa <- bsa[bsa$treat==0,]



### ideological positions (placebo)

## models

# 1984
prebsa$post <- ifelse(prebsa$year >= 1984, 1, 0)
ai_plac84 <- lm_robust(unions_lr ~ liverpool*post*lmh_l + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome,
                       clusters = constituency97_id,
                       se_type = 'stata',
                       fixed_effects = ~constituency97_id + year,
                       data = prebsa)
summary(ai_plac84)$coefficients[,1:4]

# 1985
prebsa$post <- ifelse(prebsa$year >= 1985, 1, 0)
ai_plac85 <- lm_robust(unions_lr ~ liverpool*post*lmh_l + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome,
                       clusters = constituency97_id,
                       se_type = 'stata',
                       fixed_effects = ~constituency97_id + year,
                       data = prebsa)
summary(ai_plac85)$coefficients[,1:4]

# 1986
prebsa$post <- ifelse(prebsa$year >= 1986, 1, 0)
ai_plac86 <- lm_robust(unions_lr ~ liverpool*post*lmh_l + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome,
                       clusters = constituency97_id,
                       se_type = 'stata',
                       fixed_effects = ~constituency97_id + year,
                       data = prebsa)
summary(ai_plac86)$coefficients[,1:4]

# 1987
prebsa$post <- ifelse(prebsa$year >= 1987, 1, 0)
ai_plac87 <- lm_robust(unions_lr ~ liverpool*post*lmh_l + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome,
                       clusters = constituency97_id,
                       se_type = 'stata',
                       fixed_effects = ~constituency97_id + year,
                       data = prebsa)
summary(ai_plac87)$coefficients[,1:4]


## figure

# starting
models <- list(ai_plac84=ai_plac84, ai_plac85=ai_plac85, ai_plac86=ai_plac86, ai_plac87=ai_plac87) #select models
results <- data.frame() #initialise data
variable_of_interest <- "liverpool1:post:lmh_l1" #set IV

# collect data
for (model_name in names(models)) {
  model <- models[[model_name]]
  tidy_model <- tidy(model, conf.int = TRUE)  # Calculate confidence intervals
  coefficient <- tidy_model[tidy_model$term == variable_of_interest, "estimate"]
  conf_low <- tidy_model[tidy_model$term == variable_of_interest, "conf.low"]
  conf_high <- tidy_model[tidy_model$term == variable_of_interest, "conf.high"]
  results <- rbind(results, data.frame(Model = model_name, Coefficient = coefficient, 
                                       CI_Low = conf_low, CI_High = conf_high))}

# plot
ggplot(results, aes(x = factor(Model, levels = c("ai_plac87", "ai_plac86", "ai_plac85", "ai_plac84")), 
                    y = Coefficient)) +
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = CI_Low, ymax = CI_High), width = 0.2) +
  labs(x = "", y = "Coefficient") +
  coord_flip() +
  scale_x_discrete(labels = c("ai_plac87" = "≥1987", 
                              "ai_plac86" = "≥1986", "ai_plac85" = "≥1985",
                              "ai_plac84" = "≥1984")) +
  scale_y_continuous(limits = c(-0.55, 0.55), breaks = seq(-0.4, 0.4, 0.2)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey15") + 
  theme_minimal() +
  theme(axis.text.y = element_text(size = 16), 
        axis.title.x = element_text(size = 16),
        axis.text.x = element_text(size = 14))

#save
ggsave("placebo_ideology.png", width = 5, height = 4.5, units = "in", dpi = 300)



### party support (placebo)

## models

# 1984
prebsa$post <- ifelse(prebsa$year >= 1984, 1, 0)
av_plac84 <- lm_robust(labour ~ liverpool*post*lmh_l + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome,
                       clusters = constituency97_id,
                       se_type = 'stata',
                       fixed_effects = ~constituency97_id + year,
                       data = prebsa)
summary(av_plac84)$coefficients[,1:4]

# 1985
prebsa$post <- ifelse(prebsa$year >= 1985, 1, 0)
av_plac85 <- lm_robust(labour ~ liverpool*post*lmh_l + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome,
                       clusters = constituency97_id,
                       se_type = 'stata',
                       fixed_effects = ~constituency97_id + year,
                       data = prebsa)
summary(av_plac85)$coefficients[,1:4]

# 1986
prebsa$post <- ifelse(prebsa$year >= 1986, 1, 0)
av_plac86 <- lm_robust(labour ~ liverpool*post*lmh_l + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome,
                       clusters = constituency97_id,
                       se_type = 'stata',
                       fixed_effects = ~constituency97_id + year,
                       data = prebsa)
summary(av_plac86)$coefficients[,1:4]

# 1987
prebsa$post <- ifelse(prebsa$year >= 1987, 1, 0)
av_plac87 <- lm_robust(labour ~ liverpool*post*lmh_l + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome,
                       clusters = constituency97_id,
                       se_type = 'stata',
                       fixed_effects = ~constituency97_id + year,
                       data = prebsa)
summary(av_plac87)$coefficients[,1:4]

# 1988
prebsa$post <- ifelse(prebsa$year >= 1988, 1, 0)
av_plac88 <- lm_robust(labour ~ liverpool*post*lmh_l + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome,
                       clusters = constituency97_id,
                       se_type = 'stata',
                       fixed_effects = ~constituency97_id + year,
                       data = prebsa)
summary(av_plac88)$coefficients[,1:4]


## figure

# starting
models <- list(av_plac84=av_plac84, av_plac85=av_plac85, av_plac86=av_plac86, av_plac87=av_plac87, av_plac88=av_plac88) #select models
results <- data.frame() #initialise data
variable_of_interest <- "liverpool1:post:lmh_l1" #set IV

# collect data
for (model_name in names(models)) {
  model <- models[[model_name]]
  tidy_model <- tidy(model, conf.int = TRUE)  # Calculate confidence intervals
  coefficient <- tidy_model[tidy_model$term == variable_of_interest, "estimate"]
  conf_low <- tidy_model[tidy_model$term == variable_of_interest, "conf.low"]
  conf_high <- tidy_model[tidy_model$term == variable_of_interest, "conf.high"]
  results <- rbind(results, data.frame(Model = model_name, Coefficient = coefficient, 
                                       CI_Low = conf_low, CI_High = conf_high))}

# plot
ggplot(results, aes(x = factor(Model, levels = c("av_plac88", "av_plac87", "av_plac86", "av_plac85", "av_plac84")), 
                    y = Coefficient)) +
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = CI_Low, ymax = CI_High), width = 0.2) +
  labs(x = "", y = "Coefficient") +
  coord_flip() +
  scale_x_discrete(labels = c("av_plac88" = "≥1988", "av_plac87" = "≥1987", 
                              "av_plac86" = "≥1986", "av_plac85" = "≥1985",
                              "av_plac84" = "≥1984")) +
  scale_y_continuous(limits = c(-0.55, 0.55), breaks = seq(-0.4, 0.4, 0.2)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey15") + 
  theme_minimal() +
  theme(axis.text.y = element_text(size = 16), 
        axis.title.x = element_text(size = 16),
        axis.text.x = element_text(size = 14))

#save
ggsave("placebo_support.png", width = 5, height = 4.7, units = "in", dpi = 300)





##### DD/DDD results ######

### perceived party positions

# limit to shared years
bsa_p <- bsa[!(bsa$year %in% c(1986, 1987)),]

# 1. simple DD
p1 <- lm_robust(labextr_p ~ liverpool*treat,
                clusters = constituency97_id,
                se_type = 'stata',
                data = bsa_p)
summary(p1)$coefficients[,1:4]

# 2. preferred DD
p2 <- lm_robust(labextr_p ~ liverpool*treat + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome + as.factor(lmh),
                clusters = constituency97_id,
                se_type = 'stata',
                fixed_effects = ~constituency97_id + year,
                data = bsa_p)
summary(p2)$coefficients[,1:4]

# 3. simple DDD
p3 <- lm_robust(labextr_p ~ liverpool*treat*lmh_l,
                clusters = constituency97_id,
                se_type = 'stata',
                data = bsa_p)
summary(p3)$coefficients[,1:4]

# 4. preferred DDD
p4 <- lm_robust(labextr_p ~ liverpool*treat*lmh_l + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome,
                clusters = constituency97_id,
                se_type = 'stata',
                fixed_effects = ~constituency97_id + year,
                data = bsa_p)
summary(p4)$coefficients[,1:4]

# display
texreg(list(p1,p2,p3,p4),
       caption="DD/DDD for Perceived Party Positions (Labour Extremism)",
       custom.model.names = c("1", "2", "3", "4"),
       custom.header = list("DD" = 1:2, "DDD" = 3:4),
       custom.coef.map = list("liverpool1:treat1:lmh_l1" = "Merseyside*Treat*LWC (DDD)",
                              "liverpool1:treat1" = "Merseyside*Treat (DD)",
                              "liverpool1:lmh_l1" = "Merseyside*LWC", 
                              "treat1:lmh_l1" = "Treat*LWC"),
       include.ci = FALSE,
       no.margin = TRUE,
       caption.above = TRUE,
       booktabs = TRUE,
       siunitx = TRUE)



### ideological positions

# limit to shared years
bsa_i <- bsa[!(bsa$year %in% c(1986, 1989)),]

# 1. simple DD
i1 <- lm_robust(unions_lr ~ liverpool*treat,
                    clusters = constituency97_id,
                    se_type = 'stata',
                    data = bsa_i)
summary(i1)$coefficients[,1:4]

# 2. preferred DD
i2 <- lm_robust(unions_lr ~ liverpool*treat + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome + as.factor(lmh),
                    clusters = constituency97_id,
                    se_type = 'stata',
                    fixed_effects = ~constituency97_id + year,
                    data = bsa_i)
summary(i2)$coefficients[,1:4]

# 3. simple DDD
i3 <- lm_robust(unions_lr ~ liverpool*treat*lmh_l,
                    clusters = constituency97_id,
                    se_type = 'stata',
                    data = bsa_i)
summary(i3)$coefficients[,1:4]

# 4. preferred DDD
i4 <- lm_robust(unions_lr ~ liverpool*treat*lmh_l + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome,
                    clusters = constituency97_id,
                    se_type = 'stata',
                    fixed_effects = ~constituency97_id + year,
                    data = bsa_i)
summary(i4)$coefficients[,1:4]

# display
texreg(list(i1,i2,i3,i4),
       caption="DD/DDD for Ideological Positions (Unions)",
       custom.model.names = c("1", "2", "3", "4"),
       custom.header = list("DD" = 1:2, "DDD" = 3:4),
       custom.coef.map = list("liverpool1:treat1:lmh_l1" = "Merseyside*Treat*LWC (DDD)",
                              "liverpool1:treat1" = "Merseyside*Treat (DD)",
                              "liverpool1:lmh_l1" = "Merseyside*LWC", 
                              "treat1:lmh_l1" = "Treat*LWC"),
       include.ci = FALSE,
       no.margin = TRUE,
       caption.above = TRUE,
       booktabs = TRUE,
       siunitx = TRUE)



### party support

# 1. simple DD
v1 <- lm_robust(labour ~ liverpool*treat,
                clusters = constituency97_id,
                se_type = 'stata',
                data = bsa)
summary(v1)$coefficients[,1:4]

# 2. preferred DD
v2 <- lm_robust(labour ~ liverpool*treat + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome + as.factor(lmh),
                clusters = constituency97_id,
                se_type = 'stata',
                fixed_effects = ~constituency97_id + year,
                data = bsa)
summary(v2)$coefficients[,1:4]

# 3. simple DDD
v3 <- lm_robust(labour ~ liverpool*treat*lmh_l,
                clusters = constituency97_id,
                se_type = 'stata',
                data = bsa)
summary(v3)$coefficients[,1:4]

# 4. preferred DDD
v4 <- lm_robust(labour ~ liverpool*treat*lmh_l + imp_age + female + as.factor(education) + as.factor(religion) + union + white + ownhome,
                clusters = constituency97_id,
                se_type = 'stata',
                fixed_effects = ~constituency97_id + year,
                data = bsa)
summary(v4)$coefficients[,1:4]

# display
texreg(list(v1,v2,v3,v4),
       caption="DD/DDD for Party Support",
       custom.model.names = c("1", "2", "3", "4"),
       custom.header = list("DD" = 1:2, "DDD" = 3:4),
       custom.coef.map = list("liverpool1:treat1:lmh_l1" = "Merseyside*Treat*LWC (DDD)",
                              "liverpool1:treat1" = "Merseyside*Treat (DD)",
                              "liverpool1:lmh_l1" = "Merseyside*LWC", 
                              "treat1:lmh_l1" = "Treat*LWC"),
       include.ci = FALSE,
       no.margin = TRUE,
       caption.above = TRUE,
       booktabs = TRUE,
       siunitx = TRUE)





##### Visualized ######

### prepare data

# recode
bsa$treat_plt <- dplyr::recode(as.factor(bsa$treat), "0" = "Pre-treatment", "1" = "Post-treatment")

# remove unbalanced years
bsa$labextr_p_plt <- bsa$labextr_p
bsa$unions_lr_plt <- bsa$unions_lr
bsa$labextr_p_plt[bsa$year %in% c(1986, 1987)] <- NA 
bsa$unions_lr_plt[bsa$year %in% c(1986, 1989)] <- NA

# create variable for outcomes
outcomes_df <- tidyr::pivot_longer(bsa, cols = c(labextr_p_plt, unions_lr_plt, labour),
  names_to = "variable", values_to = "value")



### calculate statistics

# initialize
results <- list()
unique_combinations <- unique(outcomes_df[, c("treat_plt", "variable")])

# calculate DDs
for (i in seq_len(nrow(unique_combinations))) {
  # extract current combination
  treat <- unique_combinations$treat_plt[i]
  var <- unique_combinations$variable[i]
  
  # subset the data for this combination
  subset_df <- outcomes_df[outcomes_df$treat_plt == treat & outcomes_df$variable == var, ]
  
  # initialize storage for subgroup means
  subgroup_means <- matrix(NA, nrow = 2, ncol = 2)
  
  # loop over `liverpool` and `lmh_l` values (0, 1)
  for (liverpool in 0:1) {
    for (lmh_l in 0:1) {
      # Subset for each subgroup
      group_df <- subset_df[subset_df$liverpool == liverpool & subset_df$lmh_l == lmh_l, ]
      
      # Calculate mean for the subgroup
      mean_value <- mean(group_df$value, na.rm = TRUE)
      subgroup_means[liverpool + 1, lmh_l + 1] <- mean_value
    }
  }
  
  # calculate DDD
  dd <- (subgroup_means[2, 2] - subgroup_means[2, 1]) - (subgroup_means[1, 2] - subgroup_means[1, 1])
  
  # store results
  results[[i]] <- data.frame(
    treat_plt = treat,
    variable = var,
    dd = dd
  )
}

# combine results
dds <- do.call(rbind, results)

# extract DDD SEs
ses <- data.frame(
  variable = c("labextr_p_plt", "unions_lr_plt", "labour"),
  ddd_se = c(
    summary(p3)$coefficients["liverpool1:treat1:lmh_l1", "Std. Error"],
    summary(i3)$coefficients["liverpool1:treat1:lmh_l1", "Std. Error"],
    summary(v3)$coefficients["liverpool1:treat1:lmh_l1", "Std. Error"]))

# merge DDs and SEs
ddd_fig <- merge(dds, ses, by = "variable")

# calculate confidence intervals
ddd_fig <- ddd_fig %>%
  mutate(pre_ci_lower = dd - 1.96 * (ddd_se/2),
         pre_ci_upper = dd + 1.96 * (ddd_se/2),
         post_ci_lower = dd - 1.96 * (ddd_se/2),
         post_ci_upper = dd + 1.96 * (ddd_se/2))



### display

# order
ddd_fig <- ddd_fig %>%
  mutate(variable = factor(variable, levels = c("unions_lr_plt", "labour", "labextr_p_plt")))

# plot
ggplot(ddd_fig, aes(x = treat_plt, y = dd, group = variable, color = variable, shape = variable)) +
  geom_errorbar(
    aes(ymin = ifelse(treat_plt == "Pre-treatment", pre_ci_lower, post_ci_lower),
        ymax = ifelse(treat_plt == "Pre-treatment", pre_ci_upper, post_ci_upper)),
    width = 0.15, size = 1, position = position_dodge(width = 0.2)) +
  geom_vline(xintercept = 1.5, linetype = "dashed", color = "black", size = 1) +
  geom_point(size = 7, position = position_dodge(width = 0.2)) +
  geom_line(size = 1.7, position = position_dodge(width = 0.2)) +
  scale_x_discrete(labels = c("Pre-treatment", "Post-treatment")) +
  scale_y_continuous(limits = c(-0.29, 0.29), breaks = c(-0.2, -0.1, 0.0, 0.1, 0.2)) +
  scale_color_manual(
    values = c("labextr_p_plt" = "#238b45", "unions_lr_plt" = "#08519c", "labour" = "#cb181d"),
    breaks = c("labextr_p_plt", "unions_lr_plt", "labour"),
    labels = c("Perceive Labour\nas Extreme", "Left-Right (Unions)", "Support for Labour")) +
  scale_shape_manual(
    values = c("labextr_p_plt" = 15, "unions_lr_plt" = 17, "labour" = 16),
    breaks = c("labextr_p_plt", "unions_lr_plt", "labour"),
    labels = c("Perceive Labour\nas Extreme", "Left-Right (Unions)", "Support for Labour")) +
  labs(y = "DDs (Residence*Class)", color = "Outcomes", shape = "Outcomes") +
  guides(
    color = guide_legend(order = 1),
    shape = guide_legend(order = 1)
  ) +
  theme_minimal() +
  theme(
    legend.position = "right",
    axis.title.x = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    axis.text.x = element_text(size = 20, color = "black"),  
    axis.text.y = element_text(size = 18),
    axis.title = element_text(size = 20),  
    legend.text = element_text(size = 18), 
    legend.title = element_text(size = 20))

ggsave("DDD_fig.png", plot = last_plot(), dpi = 300, width = 10, height = 6)


